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Nonlinear transient waves in coupled phase oscillators with inertia 
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Like the inertia of a physical body describes its tendency to resist changes of its state of motion, inertia of an 
oscillator describes its tendency to resist changes of its frequency. Here we show that finite inertia of individual 
oscillators enables nonlinear phase waves in spatially extended coupled systems. Using a discrete model of 
coupled phase oscillators with inertia, we investigate these wave phenomena numerically, complemented by a 
continuum approximation that permits the analytical description of the key features of wave propagation in 
the long-wavelength limit. The ability to exhibit traveling waves is a generic feature of systems with finite 
inertia and is independent of the details of the coupling function. 
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Coupled oscillators can exhibit a variety of dy¬ 
namical phenomena, from self-organized synchro¬ 
nization to the formation of complex phase pat¬ 
terns. In particular, coupling can enable travel¬ 
ing waves in spatially extended oscillator systems. 
Whether phase waves are possible in such a sys¬ 
tem usually depends on the details of the cou¬ 
pling function. We here show that if the individ¬ 
ual oscillators respond slowly to external coupling 
signals—a property often called ‘inertia’—, wave 
phenomena appear as a generic feature of the sys¬ 
tem, independent of the details of the coupling 
function. We investigate these wave phenomena 
numerically and analytically and show how the 
properties of the coupling function affect wave 
propagation in typical scenarios. 


I. INTRODUCTION 

Coupled oscillators are fascinating entities, not only be¬ 
cause of their robust synchronization properties but also 
as elementary constituents of pattern forming systems 1 3 . 
Spatially extended systems of coupled oscillators can ex¬ 
hibit a rich variety of patterns such as waves, solitons, 
spirals, target patterns, splay states, and other station¬ 
ary or transient phase distributions 3 11 . Among these 
patterns, traveling waves play a prominent role as they 
encode spatiotemporal information in a particularly sim¬ 
ple way. Traveling waves in coupled oscillator systems 
have been reported, e.g., in continuous reaction-diffusion 
systems 3 12 and in discrete systems, e.g., due to specific 
types of coupling 9 13 16 or spat ial gr adients of the intrin¬ 
sic frequency of the oscillators 10 17 . Depending on the 
mechanism that gives rise to wave propagation, waves 
can be persistent 81 ^ (e.g., in the case of solitons) or 
transientP^H- 2 ^ (e.g., in the case of damped waves). While 
the crucial properties that enable traveling waves have 
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mostly been sought in the particular way that oscilla¬ 
tors interact, we here consider the role of slow frequency 
adaptation of the individual oscillators for wave phenom¬ 
ena in the coupled system. This slow frequency adap¬ 
tation, often termed ‘inert ia^IIS enters the governing 
equations of a phase oscillator system as a second time 
derivative of the phase. Like the inertia of a mass point 
describes its tendency to resist changes of its velocity, 
inertia of a phase oscillator describes its tendency to re¬ 
sist changes of its dynamic frequency^. These frequency 
changes may be induced by external coupling signals or 
a time-dependent intrinsic frequency of the oscillator. 
The role of inertia for synchronization has been exten¬ 
sively studied in recent yeardP 1 31 . Inertia introduces 
qualitatively new effects such as first order phase tran¬ 
sitions in the synchronization onsetPH, cluster explosive 
synchronization^, and hysteretic transitions from inco¬ 
herence to coherence^^® 6 . We here consider the influ¬ 
ence of inertia (i.e., slow frequency adaptation) on the 
ability of a coupled oscillator system to exhibit phase 
waves. 

In this paper, we show that finite inertia enables non¬ 
linear damped phase waves in spatially extended oscil¬ 
lator systems. To this end, we study a nearest-neighbor 
coupled lattice of identical phase oscillators. This is com¬ 
plemented by a spatial continuum approximation for the 
limit of long wavelengths, which elucidates the mech¬ 
anism of wave propagation and enables to analytically 
compute the velocity and decay rate of waves as well as 
the effects of nonlinear coupling on wave propagation. 
The results show that in the presence of inertia, the abil¬ 
ity of oscillator systems to exhibit waves does not depend 
on the specific type of coupling but rather is a generic fea¬ 
ture of such systems, which enables signal propagation at 
finite velocities. 


II. COUPLED PHASE OSCILLATORS WITH INERTIA 

We study a phase oscillator model for a network of 
coupled identical oscillators with inertia, introduced by 
Tanaka, Lichtenberg, and OishP^H, based on the zero- 
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inertia model by KuramotcP, 

N 

i^Pt t Pi = si -\ ^ ^ WijV((f)j (pi ) , ( 1 ) 

J=1 

where pi is the phase of oscillator z, fi > 0 is the iner¬ 
tia of the oscillators, SI is the intrinsic frequency, k is 
the coupling strength, T is a 27r-periodic coupling func¬ 
tion, (wij) is the adjacency matrix, and rii = JA Wij 
is the total weight of links of oscillator z. We here 
consider the case of nearest-neighbor coupling on a d- 
dimensional hypercubic lattice, in which case (w ^) is di- 
agonalizable and rii = 2d. Without loss of generality, 
we consider coupling functions T that satisfy r(0) = 0 
as any finite r(0) = Tq can be absorbed by the redefi¬ 
nition £2 —>> SI + kTo and T(ip) —)> r(y?) — IV The in¬ 
ertia fi determines the time scale on which an oscillator 
adapts its frequency: for an uncoupled oscillator start¬ 
ing with initial frequency p(0) = SIq and described by 
fip + p = H, the time-dependent dynamic frequency is 
given by p(t) = SI + (SIq — Sl)e~ t ^. 

The phase waves described in the following are (pos¬ 
sibly large) perturbations of the in-phase synchronized 
solution 


pi(t) = Sit , 


( 2 ) 


in which all oscillators evolve with their intrinsic fre¬ 
quency SI and have no phase lag relative to each other. 
Hence, before turning to the description of wave phenom¬ 
ena, a comment on the stability of the solution Eq. © 
is in order. Linear stability analysis yields that this solu¬ 
tion is stable if and only if kT'(O) > 0 , a result coinciding 
with the one for homogeneous systems without inertia, 
see Appendix|Aj We only consider such cases. To simplify 
the description of waves, we consider the phase dynamics 
in the corotating frame, pi pi + Sit, which amounts to 
setting SI = 0 in Eq. ©. 

We now demonstrate that Eq. 0 can exhibit the prop¬ 
agation of damped waves for fi p 0. For simplicity, we 
here consider a one-dimensional ring of oscillators with 
nearest-neighbor coupling, = 5ij-i+Sij+i with S be¬ 
ing the Kronecker delta, and sinusoidal coupling function, 
T(0) = sin p. To show how waves propagate, we start the 
system with initial phases given by a sharply localized 
Gaussian distribution in the oscillator index z, which ap¬ 
proximately satisfies periodic boundary conditions, and 
solve Eq. 0 numerically. Fig. [l] shows the evolution of 
this initial phase distribution in a system with no inertia 
(fi = 0) and a system with finite inertia (fi p 0). For the 
case of no inertia (Fig. EH the initial phase distribution 
decays and broadens in a diffusive manner, while staying 
centered around its initial position. For the case of finite 
inertia (Fig. [l|3), the initial phase distribution initiates 
a decaying wave with finite velocity in both directions of 
the ring. In both cases, the phase differences between os¬ 
cillators eventually decay and in the limit of large times, 
the in-phase synchronized state Eq. © is reached. 



Figure 1. Phase distribution of a ring of oscillators for dif¬ 
ferent points in time for the case of (A) no inertia, fi — 0, 
and (B) finite inertia, fi — 10. Dots show a numerical solu¬ 
tion of the discrete oscillator system Eq. ([!]); curves show a 
numerical solution of the continuum approximation Eq. 0. 
Brighter curves show later times: t — 0, ...,0.5 in steps of 
At — 0.1 in A; t = 0, ..., 12.5 in steps of At — 2.5 in B. The 
space-time plots to the right show representations of the same 
systems in which each pixel along the x-axis corresponds to 
an oscillator. The color value indicates the relative phase 
of the oscillator (see color wheel in Fig. [2]4) and time pro¬ 
ceeds from top to bottom. The other parameters are N — 50, 
k = 50, and H = 0. The coupling function is r(c^) = sirup. 
Initial conditions for Eq. 0: <j>i( 0) = f{ie) and 0,(0) = 0 
and for Eq. <£(a;,0) = f(x) and d t i(a;,0) = 0, where 
f(x) — /o exp (—(x — xo) 2 /2a 2 ) with /o = 7r, xq — -£/4, and 
cr = 0.2. 


III. CONTINUUM APPROXIMATION OF THE OSCILLATOR 
LATTICE 


Wave propagation and decay can be understood us¬ 
ing a spatial continuum approximation of Eq. ([l]), which 
effectively describes long-wavelength solutions. We re¬ 
place the nearest-neighbor coupled oscillator lattice in 
Eq. ([!]) by a continuous field, pi u ...,i d (t) — > 4>(x, £), where 
(ii,..., id) G %j d and x = (eii ,..., sip) G where 6 is the 
lattice spacing. Expanding the r.h.s. of Eq. ([l]) in a power 
series in 5 to second order at e = 0 yields the nonlinear 
partial differential equation, see Appendix |B} 

(<9 2 - c 2 V 2 )<F + 2 7 <9 - A(V<F ) 2 » SI , (3) 


vhere SI = SI/fi and 




7 = V \ = f (4) 

2 Li 2d u 


Eq. |3| reveals that the Laplace operator V 2 arising 
from the expansion in the lattice spacing combines with 
the second time derivative d/ of the inertial term to a 
d’Alembert wave operator with wave velocity c. Inter¬ 
preting Eq. © as a wave equation, 7 represents a damp¬ 
ing coefficient and A describes the strength of the nonlin¬ 
ear term (V<F) 2 . For the one-dimensional ring of coupled 
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oscillators with Q = 0 considered in Fig. [T] Eq. © re¬ 
duces to 

(<9 t 2 - C 2 d 2 x )§ + 2 7 d t 3> - A(4$) 2 = o , (5) 

with periodic boundary conditions 

4>(0, t) mod 2tt = <F(^, t) mod 2tt , (6) 

d x $(0,t) = d x ${£,t) , (7) 

where is the system length. Without loss of generality, 
we here set £ = 2i r, which implies a lattice spacing 5 = 
2tt/N. We now use Eq. © to discuss wave phenomena 
of the oscillator system in the limit of long wavelengths. 


The interpretation of Eq. © as a wave equation en¬ 
ables to determine the range of the damped waves as the 
distance travelled within the time until decay as 

q 

R = — oc . (11) 

7 

Hence, the range R increases with increasing cou pling 
strength k and inertia q, while the velocity c oc yj k) q, 
Eq. © increases with k but decreases with q. In the 
limit Kfi oc with a fixed ratio ft/q, we find an infinite 
range R at finite velocity c. 

V. NONLINEAR APPROXIMATION: SPLAY STATES 


IV. LINEAR APPROXIMATION 


In Fig. EF we considered a system without inertia and 
sinusoidal coupling, implying q = 0 and A = 0. Multi¬ 
plying Eq. © by q and subsequently setting q = 0, we 
find that in this case it simplifies to a diffusion equation 

d t $ = Dd 2 x § (8) 

with constant D = c 2 /(2 7 ), which explains the 

diffusion-like decay and broadening of the initial phase 
distribution 3 ^, shown in Fig. Eh For finite inertia (q ^ 0) 
and A = 0, which corresponds to the case shown in 
Fig. [1)3, Eq. © reduces to a telegraph equation, 

£$ = ( d 2 - c 2 d 2 x + 2 7 d t )$ = 0 , (9) 


which describes the propagation of damped waves. The 
telegraph equation is linear and can be solved analyti¬ 
cally for arbitrary initial conditions 33 . We here briefly 
recapitulate the properties of damped wave solutions 
by considering the dispersion relation of a plane wave 
$(x,t) = c 1 ( kx - zt ) ^ given by Zk = \fk 2 c 2 — y 2 — iy, where 
ujk = Re Zk determines the frequency and rk = Im Zk the 
decay rate of the wave. The group velocity of a wave 
packet with wavenumbers sharply centered around the 
central wavenumber ko > y/c is given by 


n \ 

<■«(**) = -JjT 


~k 0 


( 10 ) 


which has the property v(ko) > c. Hence, such a wave 
packet propagates at least with the velocity c given by 
Eq. ©• For waves with k < y/c, Zk is purely imagi¬ 
nary and the wave is purely damped. Fig. IT shows that 
the continuum approximation (solid curves) captures well 
the propagation and decay of the waves exhibited by the 
discrete system (dots). Note that Eq. © is a linear ap¬ 
proximation of the nonlinear Eq. ([T]) in the continuum 
limit. Hence, the superposition principle that holds for 
the solutions to Eq. © only holds approximately for so¬ 
lutions of Eq. 0 with long wavelengths. 


We now turn to the full nonlinear approximation 
Eq. © with A / 0. The only non-damped traveling 
wave solutions to Eq. © of the type 

<f>(x, t) = T(x — vt) (12) 


that satisfy the boundary conditions Eqs. © and © 
splay stateJ^U, 


are 




A m 

27 


m G Z , 


(13) 


where £ = x — vt. That Eq. (13) provides the only set of 
non-damped traveling wave solutions to Eq. J5I ) can be 


seen by inserting the traveling wave ansatz Eq. ([12]) into 
Eq. ([5]), which leads to 

(v 2 — c 2 )'ip' — 2y vip — A0 2 = 0 (14) 

for the derivative 0 = 4F. This is a Riccati equation with 
constant coefficients that can be solved analytically®! 
The family of all solutions can be obtained as 

TYl 

^70 1 _ Qj^-lgmAy^-c 2 ) 5 (1^) 

where m = — 2jv/\ and a E [— 00 , 00 ]. The boundary 
condition Eq. 0 translates into 0(£) = 0(£ + £). How¬ 
ever, the only ^-periodic members of the family 0 a are 
the constant solutions 0o = m and 0oo = 0. The bound¬ 
ary condition Eq. © forces m to be integer. Hence, the 
only family of solutions for T is given by Eq. (13). Since 


the second time derivative of the splay states dynami¬ 
cally vanishes, 4> = v 2 T^ = 0, they coincide with the 
splay states found in systems without inertia (q = 0). 

To show that Eq. © describes the long-wavelength 
splay states of the discrete model, we use the splay state 
ansatz <fii(t) = m(si — vt) in Eq. (|T]) , which leads to 

2 mv _ N 

—— = T(em) +r{-em) 

K (16) 

— £ 2 m 2 T"{0) + 0{e A m A ) . 

This yields v ~ K£ 2 mr"( 0)/2, in accordance with 


Eqs. (13) and ©. Note that this result for the con¬ 
tinuum approximation Eq. © does not rule out the ex¬ 
istence of other non-damped traveling wave solutions to 
Eq. 0 that may arise from higher order nonlinearities 
in the coupling function T. 
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VI. NONLINEAR APPROXIMATION: DAMPED WAVES 

We now investigate how the effects of the nonlinearity 
in Eq. © modify the propagation of damped waves. We 
here consider cases where the nonlinearity yields a small 
perturbation to the solution of the linear telegraph equa¬ 
tion © and use a perturbative expansion in A. Using 
the ansatz 


4> = 4 > 0 + A4>i (17) 

in Eq. (J5|, we obtain, to zeroth and first order in A, the 
two linear equations 


£4 > 0 = 0 , (18) 

£$! = (d x $o) 2 , (19) 


where the linear operator C of the telegraph equation was 
defined in Eq. ([9|. In the following, we use Eqs. (17-19) 
to study two generic cases for wave propagation: delocal¬ 
ized waves extending throughout the system and a single 
localized wave. 


A. Delocalized waves 


As a first generic scenario we consider the behavior 
of delocalized waves that extend throughout the system. 
As a starting point for the perturbative treatment, we 
consider the following exact solution to Eq. (18) with 
boundary conditions © and 0 , 


4>o(x, t) = 0 cos {kx — lO kt)e 


- 7 1 


( 20 ) 


where fc E Z, uJk = \/k 2 c 2 — y 2 , and (j) is the amplitude 
of the wave. Eq. (20) describes a damped plane wave 


with wavenumber &, phase velocity Vk = c Ok/k and decay 
rate 7 . We only consider cases in which is real, that is, 
kc > 7 . Eq. (19) can be solved analytically using Fourier 
transformation, see Appendix [Cj 

A$i(M)= [l-(l + 2 7 t)e- 2 ^+p(x,t)e-^]A , (21) 
where 


A = 


Xcj) 2 k 2 
87 2 


( 22 ) 


The term p(x, t) contains higher harmonics in the 
wavenumber k but remains bounded as 


\p(x,t)\ <— + V2e 


(23) 


for all cases in which the frequency ujk is real, see Ap¬ 
pendix [Cj Since Eq. ( 21 ) implies <f>i —» A as t 00 , the 
asymptotic effect of the nonlinearity is a global phase 
offset of magnitude A. Hence, for delocalized waves, the 
nonlinearity effectively acts as a global source (or sink, 


depending on the sign of A) that transiently contributes 
to the global frequency. 

To demonstrate that our analytical results present an 
effective description of the discrete system Eq. 0 in the 
limit of long wavelengths, we compare them to numeri¬ 
cal solutions of Eq. 0 for different wavenumbers k and 
different strengths A of the nonlinearity. For the discrete 
model, we choose different coupling functions from the 
family 


r A (v?) = sin<£ - A(cos^ - 1) 


(24) 


where we have identified the parameter A of the contin¬ 
uum approximation according to Eqs. 0 , see Fig. §\. 
Fig. [2p> shows a comparison of solutions to the discrete 
system Eq. ( [Th (u pper row) and the analytical approxi¬ 
mation, Eqs. (|17|), (20) and ( | 2 l| ) (lower row), for differ¬ 
ent values of A. Compared to the analytical approxima¬ 
tion, the discrete model shows undulations in the width 
of wave peaks which are the result of higher order non- 
linearities of the coupling function (first three panels in 
Fig. [2^3). For weak nonlinearity (small A), the analytical 
approximation captures the time evolution of the system 
very well (first three panels in Fig. [2^3), while for stronger 
nonlinearity (large A), the approximation breaks down 
(rightmost panel in Fig. [2^3) as the discrete system ex¬ 
hibits more complex phase patterns. Fig. [2p shows the 
same comparison for fixed A but different wavenumbers 
k. For longer initial wavelengths (small fc), the results 
of the discrete model and the continuum approximation 
agree better and better (first three panels in Fig. [2p). 
For small wavelengths (large fc), the long-wavelength con¬ 
tinuum approximation breaks down (rightmost panel in 
Fig. [2p). In all cases considered here, the system even¬ 
tually reaches the in-phase synchronized state Eq. ©• 


B. Localized waves 


As a second generic case we consider the propagation 
of a single localized wave on an infinite spatial domain. 
Hence, we drop the boundary conditions Eqs. © and 
0. As starting point for the perturbative treatment 
Eqs. (17-19), we use a traveling Gaussian wave packet 


that decays, 


$0 (x,t) = 4>e 


= rhp -( x ~ ct ') 2 / 2ff2 e -7 * 


(25) 


initially centered around x = 0 with width < 7 , height 0 , 
velocity c, and decay rate 7 , see Fig. (red curves). 
Eq. (25) is an approximate solution to Eq. (H8l) for weak 
damping, £<f >0 = 0( 7 2 ). Using a Fourier transformation 


in the spatial domain, an exact solution to Eq. (19) can 


be obtained for 4>i (k,t) = / <&i(x, t) e~ lkx dx. However, 
the corresponding backtransform to real space cannot be 
represented in terms of elementary functions. Hence, we 
approximate the exact solution for 4>i (fc,t) in terms of 
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Figure 2 . Time evolution of delocalized waves in a ring of oscillators. (A) Coupling functions Ta, Eq. ( p4| , for A = 0 
(solid), A = —0.2 (dashed), and A = 0.2 (dotted). (B) Numerical solutions of the discrete model Eq. 0 (upper row) and the 
approximate continuum solution (lower row), Eqs. fl7| , ( [20| and (21), where the contribution of p(x,t) has been neglected, 
for different A but same wavenumbers k. (C) The same as in B but for fixed A and different wavenumbers k. In B and C, 
the coupling function of the discrete model is given by Eq. (24). The other parameters are p = 6 , k = 20, = 0 for Eq. 0 , 

implying c ~ 0.16 and 7 ~ 0.08 for Eq. 0. Initial conditions for Eq. 0: rf>, (t ) = < 1 >q('<c, 0) and <i,(t) = (VlVifis. 0). and for 
Eq. |5|: $(x,t) — 4>o(x,0) and dt$(x,t) = <9t4>o(:r, 0), with 4>o given by Eq. (|20|) with 0 = 37t/4. 


Gaussian and trigonometric functions, which permit a 
backtransformation, see Appendix |Dj We thus obtain 


(x, t) ~ (j) 2 [(1 — e ' yt )U(x — ct) + T (x, t ) 
+ \h(x + ct) — \h(x — c£)]e _7t , 


(26) 


where the functions II, \h, and T are given by 


U{X) = ~A^ e{x,a) ’ 

*(x) = I^[ e_Wcr)2 +2(x/a)6(x/a)} , 

Y(x,t) = 1 4 -yuc 


(27) 

(28) 

(29) 


with 


6(x) = f e y dy 

Jo 

being an error function and 


= 


i 

\J\ + 2c 2 t/7<r 2 


(30) 


(31) 


being a time-dependent scaling factor. Eq. ( |26j ) reveals 
that the effect of the nonlinearity can be understood by 
three contributions: since 6 essentially is a smoothed step 
function, the contribution II (x — ct) represents a front 
with step width <j that propagates with the center of the 
wave. The prefactor 1 — e _7t describes a gradual build¬ 
up of this front. The behavior of the contributions 4/ 
and T is more complicated but effectively deforms the 
front into a tail that the traveling wave leaves behind; 
this behavior of <f>i and the full perturbative approxi¬ 
mation Eq. © is displayed in Fig. [3]A (green and blue 
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A = 0, <t = 0.4 A = 4 x 1CT 3 , o- = 0.4 A = -8 x 10" 3 , <r = 0.4 A = -1 x 1CT 2 , <r = 0.4 



1 / 50 




Figure 3. Time evolution of a localized wave on an infinite domain. (A) Localized decaying wave <t>o (red), Eq. (25), the 
perturbative contribution A<f>i (green), Eq. (26), and the sum (blue), Eq. ( |17| ). Brighter curves correspond to later times 
(i — 0,..., 60 in steps of At = 4). Parameters as in the second panel from the left of B. (B) Numerical solutions of the discrete 
model Eq. ([Tj) (upper row) and the approximate continuum solution (lower row), Eqs. |l7| , (25) and (26), for different A but 
same width a of the wave. (C) The same as in B but for fixed A and different widths a of the wave. The other parameters are 
/i = 10, k = 50, Q = 0 for Eq. 0 , implying c ~ 0.2 and 7 = 0.05 for Eq. 0 . Initial conditions analogous to Fig. |2^,C but 
with $0 given by Eq. (25) with <j) — 37 r/ 4 . 


curves). Figs. [3|3,C compare numerical simulations of 
the discrete model Eq. 0 to this analytical approxima¬ 
tion for different A and different widths a of the trav¬ 
eling wave. Again, for weak nonlinearity (small A) and 
long wavelengths (large cr), the analytical result ( [26] ) ap¬ 
proximates the discrete model very well (first three pan¬ 
els in Figs. [3^3, C). For strong nonlinearity (large A) and 
short wavelengths (small cr), the continuum approxima¬ 
tion breaks down and the higher order nonlinearities in 
the coupling cause the tail of the wave to evolve into com¬ 
plex phase patterns (rightmost panels in Figs. [3j3,C). 


VII. DISCUSSION 

In this work, we have shown that the interplay of in¬ 
ertia and coupling enables traveling waves in spatially 
extended oscillator systems, independent of the details 
of the coupling function. A spatial continuum approxi¬ 


mation for long wavelengths has been derived that elu¬ 
cidates the mechanism of wave propagation and permits 
the analytical computation of key observables such as the 
wave speed and damping rate. We found that the range 
and velocity of these waves is governed by the magnitude 
of inertia and coupling strength. For the case of purely 
sinusoidal coupling, wave propagation can be described 
by a linear telegraph equation to a very good approxi¬ 
mation. For general coupling functions, the continuum 
approximation is nonlinear. This nonlinearity modifies 
the behavior of damped waves and enables splay states 
with finite frequencies corresponding to persistent trav¬ 
eling waves. For delocalized waves, we found that the 
nonlinearity transiently modifies the global frequency of 
the system; for localized waves, the nonlinearity leads to 
a tail that the traveling wave packet transiently leaves 
behind. 

The continuum approximation Eq. 0 is a generaliza¬ 
tion of many well-known equations. It generalizes the 
telegraph equation 0 describing damped linear waves 
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to include nonlinear effects. It is a generalization of the 
Kardar-Parisi-Zhang (KPZ) equation 36 without noise to 
include a second time derivative. Differentiating Eq. © 
with respect to x and defining u = d x & yields a general¬ 
ized Burgers’ equation d 2 u + 2y d t u = c 2 d 2 u + A ud x u, 
which contains an inertial term. For A / 0 and ap¬ 
propriate boundary conditions, this generalized Burgers’ 
equation is known to exhibit non-decaying traveling wave 
solutions 37 . Moreover, equations of the type Eq. (J3J) have 
been used to describe the propagation of electromechan¬ 
ical waves in large electric power system^I 39 l 

The results presented here show that the interplay of 
coupling and inertia (i.e., slow frequency adaption) en¬ 
ables signal transmission and the distribution of spatial 
information through phase waves. These findings may 
therefore be relevant for natural and engineered cou¬ 
pled oscillator systems with slow frequency adaptation in 
which phase patterns can play a functional or interfering 
role, such as in biological oscillators^, power grids 40 4 ®, 
and systems of coupled electronic clocks 42 . More gener¬ 
ally, the present work highlights how the internal dynam¬ 
ics of an oscillator can crucially contribute to the emer¬ 
gence of collective phenomena that appear in coupled 
systems. While we here investigated phase waves with 
long wavelengths, more general phase dynamics will un¬ 
fold if the long-wavelength criterion is abandoned. This 
remains a rich subject for future research. 


where Wij = Wij/rii is the normalized adjacency matrix. 
To decouple this set of equations, we define the collective 
modes $i = ^ .(ra -1 )^^-, where (m^-) is the matrix that 
diagonalizes (wij) according to ^Z k i(m~ 1 )ikWkirrikj = 
diag(ei,..., ejv) and are the eigenvalues of w. This 
yields the decoupled set of equations /idi + gi$i = 0, 

where gi = kT'(0)(1 — e$). The exponential ansatz $i(t) = 
e Zit yields the characteristic equation fiz 2 + Zi + gi = 0 
with solutions 


*»,± 


—l =b yT — 4 figj 
2 g 


(A2) 


The in-phase synchronized state is linearly stable if and 
only if Re z^± < 0 for all i. Using Gershgorin’s circle the¬ 
orem, it has been shown 43 ^ that the eigenvalues of the 
normalized adjacency matrix wij satisfy |e*| < 1. This 
implies sign gi = signKT'(O). Consequently, R ezi : ± < 0 
for all modes i if and only if kT'(0) > 0, in which case 
the synchronized state is linearly stable. This is the same 
criterion as for systems without inertia 43 . 
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Appendix A: Linear stability analysis of the in-phase 
synchronized state 

To assess under which circumstances the in-phase syn¬ 
chronized state Eq. © is stable, we use the ansatz 
(j)i(t) = fit + g(i(t) in Eq. ([l]) , where Q is a perturba¬ 
tion of order unity and g is a small expansion parameter. 
Linearizing the resulting dynamics in g yields the gov¬ 
erning equations for the perturbations as 

N 

Ki + Q = «r'(o) Y, *>a (0 - C i), (Ai) 

0 = 1 


Appendix B: Derivation of the spatial continuum 
approximation Eq. (3) 

In this section, we derive the spatial continuum approx¬ 
imation Eq. © from the discrete model Eq. 0 with 
nearest-neighbor coupling. We proceed along the lines of 
Refs.® 32 . We first show the derivation for d = 1 dimen¬ 
sion and then extend the result to arbitrary dimension¬ 
ality d. The one-dimensional nearest-neighbor adjacency 
matrix is given by = Sij- i + We replace (j)i(t) 

by <£>(£,£) and 4>i± a (t) by =b crs,t) with a E { — 1,1} 
in Eq. (1), where 6 is the lattice spacing, 


+ d t <s> = n + £ J2 r ($(* + ae ’ *) - $(*> *))» 

<T = ±1 

(Bl) 

and expand the coupling function in the lattice spacing 


2 2 

r ($(x + ae,t) - $(x,t)) = a£r'(0)d x <f>(x,t) + i-r'(0)flg$(s,f) + jT"(0)[d x $(x,t)} 2 +0(£ 3 ) , (B2) 











where we have used a 2 = 1. Using the expansion Eq. (B2) in E q. (Bl), we find that terms of linear order in 6 cancel 
as the cr-sum runs over the values 1 and —1. Dividing Eq. (Bl) by fi and dropping terms of order e 3 , we obtain 


# 2jdt& = (i + c 2 <9 2 $ + A (d x $y , 


(B3) 


where Q = Q//i, and 7 , c, and A have been defined in Eqs. (|4j). The analogue of Eq. (Bl) for arbitrary dimension d is 

d 

+ d t * = fi + ^ E E r($(x + aee n ,t)-$(x,t)) , (B4) 

n= 1 cr=±l 

where e n is the unit vector in n-direction. The corresponding expansion of the coupling function is given by 

f )<T> <r 2 p 2 / \ 2 

T($(x + cr£:e n ,£) - $(x,*)) = aer'(O) —(x,J) + y r '(0)^(x,t) + yT"(0) ^ —(x,t)J + 0(£ 3 ) . (B5) 


Dividing Eq. (B4) by fi and using the expansion Eq. (B5) 
yields 


d r 


2 <9 2 $ ( d$ 

dxl + V<9x n 


(B 6 ) 


d 2< h + 27<9 t <h = D + ^ 

n=l L n 

= U + c 2 V 2 4> +A(V4 >) 2 . 
This completes the derivation of Eq. ©• 

Appendix C: Derivation of Eqs. {21} and {23) 


In this section, we derive Eq. ( |2T} , that is, we solve 
Eq. (19) with 4>o given by Eq. (1201) . For convenience, 


we here set the amplitude 0=1 and restore it later. The 
rhs of Eq. (19), given by (d x <&o) 2 = k 2 sin (kx—uJkt) 2 e~ 2lt 
with Uk — \/k 2 (? — 7 2 , can be rewritten using standard 
trigonometric relations, 


(4^o ) 2 = y 


1 — cos( 2 kx) cos( 2 Wkt) 
— sin( 2 kx) sin( 2 Ukt) 


(Cl) 


^-2 it 


I 

Hence, we make the ansatz <hi(x,0 = ipo(t) + 

cos(2kx)ipi(t) + sin( 2 /cx)(^ 2 (^) in Eq. ([19} , which yields 
the following ordinary differential equations for <^q, <£ 1 , 
and </? 2 , 


<Po + 27(^0 = ye 274 , 
fc 2 

01 + 4ft 2 c 2 ^ + 27^1 = cos( 2 u> k t)e~ 2ryt , (C 2 ) 

k 2 

02 + 4ft 2 cV 2 + 2702 = - y sin(2u>*0)e 2jt . 


Solving Eqs. ( |C 2 | ) with the initial conditions (/?o(0) = 
^i(O) = ^ 2 ( 0 ) = 0 leads to Eq. ( 21 ), where p(x,t) is 
given by 


-lp(a;,f) = (sinfh^t) + — cos (uj 2 kt) ) sin( 2 kx) - 
A k \U2ki 7 


2k 2 c 2 -Y . 


7^2/c 


sin(cj 2 /c^) — cos^/cO cos (2 kx) 


sin( 2 Ukt) + — cos( 2 uo^t) ) sin( 2 kx) — ( — sin( 2 u;fct) — cos( 2 uo^t) ) cos( 2 kx) 


^- 7 1 


(C3) 


with Afc = 7 /fee. We now prove the bound (23). Us¬ 
ing the triangle inequality \a + b\ < \a\ + \b\ and the 


trigonometric inequality \a cosx + 6 sinx| < V a 2 + 6 2 , we 
successively obtain the bounds 
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p-|pOM)l < 


Uk 


Uk 


sin(cj 2 /c^) H-cos(u; 2 /c^) ) sin( 2 kx) — 

^2k 7 


2k 2 c 2 — 7 2 
7^2/c 


sin(u; 2 /c^) — cos(cj 2 /c^) ) cos( 2 kx 


< 


< 


sin( 2 cj/ c t) + — cos(2uy7) ) sin( 2 &x) — ( — sin( 2 co’/ c t) — cos( 2 uo^t) ) cos( 2 kx) 
7 ) V 7 / 

Wk . , , cj/c , A 2 (2k 2 c 2 — 'y 2 . \ 2 i i / 2 

-sin(w2fet H- cos(oj 2 kt) + -sm(cj 2 /c^) - cosm^G 

^2/c 7 /V 7 ^2 /c 


3 - 7 t 


+ 


sin( 2 u;fct) + — cos( 2 ujkt) ) + ( — sin( 2 cc;^t) — cos( 2 uj^) 
7 ) V 7 


2 -i !/ 2 




, / 2 A: 2 c 2 — 7 2N 2 
7 


w 2 ik / v 7 y v 7^2* 
/ l — A|/2 
IAfc | \ V 1 Afe/4 


+ e 


-7t 


1/2 


2 + 21 — 


1/2 


a ~7* 


Hence 


! , |p(a:,t)| < y/2b(A k ) + \/ 2 |A fc |e 7t , where 


6 (A) = | Ah 


7-A 2 /2 
1-A 2 /4 ' 


(C4) 


where <//,■ = fee/ 7 and the functions W a , f/ a: . and 1 + are 
given by 

Wa(»)= Q-“V)e"“V (D3) 


We only consider cases in which the w ave is not purely 
damped, which corresponds to cj/c = y/ k 2 c 2 — y 2 being 
real. This implies |fcc| > |y| or, equivalently, |A&| < 1. 
The bound (23) then follows from the observation that 


0 < 6 (A) < V2/3 for |A| < 1 . 


and 


Ua(y) = 


Va(y ) = 


y.y/l-y 2 _ -a^/l^y 


2a,\J\ - y 2 


Ay/l-y 2 i -ayA 


(D4) 

(D5) 


Appendix D: Derivation of Eq. 


In this section, we derive the approximate solution 
Eq. (26). We start from Eq. ( fl9| ) with <l>o given by 
Eq. (25). As in Appendix [Cj we set the amplitude 
4> = 1 for convenience and restore it later. Express¬ 
ing <I>i by its spatial Fourier transform as <hi (x,t) = 
(2?r)~ 1 f e lkx &i(k,t) dk, Eq. (19) can be written as 


9 2 4>i + & 2 c 2 4>i + 2yd t <&i = J dxe 
1 k 2 a 2 


—ikx 


(^$o ) 2 


V?r 
o V 2 


e -fe 2 <r 2 /4 e -(2 7 +ifcc)t 


(Dl) 


Eq. ( |D1| ) is a linear ordinary differential equation in t 
whose exact solution is given by 


$i(M) = 


2y 2 cr 


(yk i')Uyt(yk) r Yt T t{yk) 

- ie - ( iyfc + 1 ) 7 A ^7 cr/2c(l/k) 

/ Vk 


(D2) 


a . Approximation of W a The function W a 


Eq. (D3), can be approximated by ~ W a , where 


Wa(y> = ^e-“ 2 ^ 2 , 


(D 6 ) 


since the term — a 2 y 2 becomes of the same magnitude 
as the term 1/2 at y = l/y/2a\ this is of the same 
order as the decay length 1/a of the Gaussian term 
e~ a y , which thus suppresses this contribution. The 
maximum distance between W a and W a is given by 
D w = sup y | W a (y) — W a (y)\ = whic h de - 

creases with decreasing a. Since a oc a in Eq. (D 2 ), 


the approximation becomes more accurate for sharper 
localized phase perturbations. 

b. Approximation of U a and V a Due to the expres¬ 
sion v^- y 2 in the exponentials, the functions U a and 
V a do not possess a Fourier transform in terms of ele¬ 
mentary functions. Hence, we approximate U a and V a 
by Gaussian and trigonometric functions, which permit 
a backtransform of 4>i to real space. We construct ap¬ 
proximating f unction s in the following way. First, note 
that, while y/l — y 2 may become imaginary for \y\ > 1 , 
both U a and V a are real. For \y\ < 1, this is obvious, 
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while for \y\ > 1, we can rewrite U a and V a as 

e i a Vv 2 ~i — e _1 u^/y 2 - 1 sin^y ^y 2 — 1) 


Ua(y) 

Va(y) 


12/1 > 1 


12/1 > 1 


2ia v /y 2 - 1 

Jay/y^l , g-iay'y 2 -! 


a 




cos(a^/y 2 — 1) • 
(D7) 


For large \y\, we approximate ^/j/ 2 — 1 ~ |y| and thus 
obtain the asymptotic behavior 


tfa(y) 

K,(y) 


ls/l>l Oi\y\ 

~ cos(ay) . 

hl>i 


sin ay 
ay 


(D8) 


Both f7 a and W exhibit a global maximum at y = 0, see 
Fig. [4|4,B. For small |y|, we approximate this feature of 
U a and V a by Gaussian functions. To obtain approxi- 
mati ons c onsistent with the asympotic behavior given by 
Eqs. (D8), we make the ansatz 


U a (y) = e u ^ + S ^ 


ay 


V a (y) = e Va(y ) + cosay , 


(D9) 


where the functions u a and v a are to be determined in 
such a way that they capture the behavior of U a and 
V a for small \y\. T o thi s end, we solve Eqs. (D9) for u a 
and v a and, using (D4) and (D5), expand the resulting 
expressions in y. 


Ay) = log U a (y) - 


sin ay 
ay 

/sinha \ p a 2 3 

= ° S \ a - 1 ) ~ T V + °( y ’ 


v a (y) = log (V a (y) - cos ay) 

= log (cosh a - 1) - y y 2 + 0(y 3 ) 


(DIO) 


(Dll) 


where 


Pa = 


a + a 3 /3 — a cosh c 
a — sinh a 
1 a sinh a — a 2 


- 1 


(D12) 


Qa = a 


2 sinh(<a/2) 2 


Hence, we obtain the approximations U a ^ U* and V a ~ 
( sinh<a 


V*, where 


Kiy) = 


\ a 


1 e -p«v / 2 + 


sin ay 
ay 


(D13) 


V*(y) = (cosha - l)e 


_ i \p-Qc*y 2 /2 


+ cos ay . 


Note that p a and q a are complicated functions of a which 
lead to a com plica ted time dependence of the approxima¬ 
tion for Eq. (D2), where a = jt. We can considerably 


simplify the approximations for U a and V a by further 
approximating U^ — U a and V* ~ V a , where 


/sinha - \ ay 2 / 2 smay 


U a (y) = I-lle-“"^ + - , / n 

V a ) a V (D14) 

V a (y) = (cosh a — l)e~ ay + cos ay . 


Justifying the approximations (D14) is not straightfor¬ 
ward as they do not solely rely on an approximation of 
p a and q a . Hence, we first give some intuitive reasoning 
and afterwards show that the approximations are valid 
by numerically com puting the relative distance bet ween 
the appro ximations (D14) and t he ex act functions (D4) 
and ( |D5 
observations: 


The approximations (|D14j) are based on two 
(i) for a 


oo, p a and q a asymptotically 
behave as p a ~ a — 1 and q a ~ a, which is straighfor- 
ward to show, and (ii) the prefactors (<a -1 sinha — 1) = 
a 2 /6 + 0(a 4 ) and (cosh a — 1) = a 2 /2 + 0(a 4 ) 


m 


Eqs. ( |D14 ) suppress the effects of the Gaussian terms 
for small a and eventually vanish for a —> 0. Hence, 
our approximation strategy is to replace p a and q a by 
their asymptotic behavior at large a. The effects of the 
large error for small a that is thus introduced is then sup¬ 
pressed by the prefactors. While the replacement of q a by 
a leads to no difficulties at this stage, the replacement of 
p a by a — 1 would lead to positive exponents in the Gaus¬ 
sian for a < 1, which would consequently lead 

to an unbounded growth as \y\ —>• oo. Hence, we approx¬ 
imate a — 1 ~ a, which is a reasonable approximation for 
a 1, while still leading to a bounded behavior of U a for 
a < 1. Again, the effects of the large error introduced for 
small a remains suppressed by the prefactor. Note that 
the thus introduced deviations from the exact values p a 
and q a only alter the width of Gaussians and not their 
height; for instance, the exact identities U a ( 0) = U a ( 0) 
and V a (0) = V^(0) are preserved by this approximation. 
Fig. [4|^,B show co mpar isons of the exact function U a and 
V a , Eqs. ( [D4| a nd (D5), with the approximations U a and 
V a , Eqs. (|D14|), for different values of a. To show that 


that Eqs. (D14) provide reasonable approximations, we 


numerically compute the relative distances 


T^u (a) = 
IV (a) = 


\\U a -U a \\ 


OO 


\\U< 

\\V a ~ V c 


a ||oo 
a ||oo 


(D15) 


\\Vc 


a oo 


where ||/||oo = sup^ \f(y)\ is the uniform norm. The 
result is s hown in Fig. [4p highlighting that the approx¬ 
imations (D14) work especially well for very small and 


very large a, while satisfying Vjj < 15% and Vy < 15% 
for all a in the computed range. In summary, the ap- 
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A 





Figure 4. (A,B) Compari son of the functions U a and V a (solid curves), Eqs. (D4) and (D5), with the approximations U a and 
V a (dashed curves), Eqs. (D14), for a = 0.5 (green), a = 1 (gray), and a = 2 (black). (C) Relative distances T>u (yellow) and 


^es) _ 

TV (green), Eqs. ( |D15| ), as a function of a. 


proximation of 4>i is given by 


$l(M) =* 


27 2 <j 


(yk i)Uyt(yk)ift t 

) 


— ie - ( iy/c + 1 ) 7t ^ W'yg/2c(yk) c _ 7 t 
Vk 


(D16) 


where W a , U a , and V a are given by Eqs. (D6) and 
(D14). The corresponding Fourier transform to real 
space, <Fi (x,t) = (27t) _1 f e lkx &i(k, t) dfc, is given by 

Eq. @. 
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